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INTRODUCTION 


Sets  of  point-masses  have  been  recognized  in  recent  years  as 
viable  models  for  describing  the  anomalous  gravity  field  of  the  earth.1'2 
The  value  of  point-masses  stems  from  the  fact  that,  in  particular 
circumstances,  computation  of  the  anomalous  potential,  its  gradient  (the 
gravity  disturbance  vector) , or  other  related  quantities  are  more 
rapidly  achieved  by  using  point-masses  instead  of  various  integral 
equation  methods  or  spherical  harmonic  expansions.  Point-masses  have 
also  proved  useful  in  describing  local  variations  of  planetary3  '**  and 
lunar5  gravity  fields. 

In  this  report,  emphasis  will  be  placed  on  the  use  of  point-mass 
sets  to  describe  the  global  anomalous  gravity  field.  First,  the  spher- 
ical harmonic  expansion  of  a point-mass  set  will  be  developed.  This 
will  allow  the  definition  of  "global  errors"  (in  the  potential  and 
in  its  gradient)  which  arise  when  a finite  spherical  harmonic  expansion 
is  replaced  by  a set  of  point-masses.  Then,  exact  forms  of  these  global 
errors  will  be  determined. 

The  development  of  analytic  expressions  for  the  global  error  in  the 
anomalous  potential  and  its  gradient  will  permit  the  establishment  of  a 
means  by  which  a point-mass  set  may  be  optimized.  Defining  an  "optimal" 
point-mass  set  as  one  which  minimizes  global  gradient  error  will  lead  to 
a system  of  optimizing  equations.  The  structure  of  these  equations,  as 
well  as  methods  for  their  solution,  will  also  be  discussed. 

Thus,  the  primary  concern  of  this  report  is  to  provide  a theoret- 
ical basis  for  the  construction  of  optimal  point-mass  sets  from  a 
finite  spherical  harmonic  expansion.  Point-mass  sets  can  be  construct- 
ed, of  course,  by  methods  which  do  not  explicitly  depend  on  spherical 
harmonic  representations;  however,  since  the  coefficients  of  a spherical 


1 


harmonic  expansion  can  be  derived  from  many  sources,  including  satellite 
orbits,  gravity  anomalies,  and  geoid  heights,  the  results  presented  here 
provide  a general  prescription  for  using  a knowledge  of  any  of  these 
basic  sources  to  create  an  optimal  point-mass  set. 

THE  SPHERICAL  HARMONIC  EXPANSION 
OF  A POINT-MASS  SET 

An  important  purpose  of  point-mass  modeling  is  to  provide  a de- 
scription of  the  variations  in  gravity  close  to  the  surface  of  the 
earth.  These  variations  are,  in  general,  due  to  mass  distributions  both 
inside  and  outside  the  earth;  however,  since  external  masses  are  easily 
identified  (the  atmosphere,  the  moon,  etc.),  their  effects  can  readily 
be  dealt  with  (i.e.,  subtracted  from  geophysical  measurements),  leaving 
data  which  reflects  only  the  influence  of  internal  masses.  (In  this 
report,  the  term  "geophysical  data"  will  be  used  to  represent  these 
subtracted  measurements.)  Thus,  the  point-mass  sets  under  consideration 
will  all  be  internal  to  the  earth. 

The  anomalous  potential  is  usually  defined  as  the  small  difference 
between  the  actual  potential  and  some  suitably  chosen  reference  po- 
tential.6 In  the  ubiquitous  "spherical  approximation"*  it  may  be 
represented  by  a spherical  harmonic  expansion,  an  expansion  whose  coef- 
ficients can  be  derived  from  a knowledge  of  gravity  anomalies,  geoid 
heights,  or  other  sources.  The  anomalous  potential  T may  also  be 

!! 


* See  Reference  6,  p.  87 
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represented  by  a set  of  point-masses,  in  which  case  it  takes  the  form 


T (r) 


I 


-k  E 

i=l 


(1 


where 


k:  the  gravitational  constant 

: the  value  of  the  ith  point-mass 
?i:  the  position  vector  of  the  ith  point-mass 
r:  the  position  vector  of  the  computation  point 
I:  total  number  of  point-masses 


The  geometric  relation  between  r and  one  of  the  ri  is  shown  in 
Figure  1.  In  that  figure,  P is  the  computation  point  and  P^  is  the 
location  of  the  ith  point-mass. 


Figure  1.  Geometry  of  Position  Vectors 


The  coordinates  of  the  points  P and  are  (r,  <j>,  X)  and 

(rj.,  <(>i/  Xi) , i -spectively,  in  a geocentric  coordinate  system  where  r 
is  the  radial  distance  from  the  center,  4>  is  the  latitude,  and  X is 
the  longitude.  In  Figure  1,  i|<i  is  the  angle  between  the  vectors  r and 
r^  and  the  relationship  between  and  the  geocentric  coordinates  of 
r and  r^  is 

cos  = sin<[>  sin<t>i  + cos<)>  cosi{>i  cos(X-Xi)  (2) 

The  distance  between  the  points  P and  Pi  is  £i,  where 

Jli  = |r-ri|  = (r2-2rricos^i  + r^)*5  (3) 

The  potential  T can  be  written  in  a spherical  harmonic  expansion 
by  using  (1)  along  with  the  following  results.  First,  the  expansion  of 
in  Legendre  polynomials 


4r  = £ 


oo  n 

ri 


Pn (cos  ^i) 


Second,  the  expansion  of  the  Legendre  polynomials  Pn(cos  i|>)  in  terms  of 
geocentric  coordinates  (the  "addition  theorem") 

n 

Pn  (cos  i^i)  = Y,  Nnm  pnm<sin  <t>)Pnm(sin  ^i^  * 

m_°  (cos  mX  cos  mXi  + sin  mX  sin  mXi) 

where  the  Pnm  are  associated  Legendre  functions  and 


l (n~m)  ! 

Nnm  em~  (n+m)  ! 


[ 1,  m=o) 


Now,  if  (5)  is  put  into  (4)  and  the  result  of  this  is  put  into  (1) , 
then  the  following  spherical  harmonic  expansion  ensues: 


n=o  m=o 


n+1 


T(rj  =~^  H £ (f)  pnm  (sin  <f>)  |Anm  cos  mX  + Bnm  sin  ml 


(7) 


where 


'nra 


v-  Mi 


pnm 


(sin  4>i) 


I cos  mX^j 
sin  mXjJ 


(8) 


(Here,  R:  mean  earth  radius,  and  M:  earth  mass.  These  parameters  ap- 
pear in  (7)  and  (8)  in  order  to  make  the  coefficients  Annl  and  Bnm 
dimensionless. ) 


Expression  (8)  shows  explicitly  how  to  construct  the  spherical 
harmonic  expansion  equivalent  to  a given  point-mass  set.  The  main 
concern  of  point-mass  modeling,  however,  is  the  converse  problem: 
given  a spherical  harmonic  expansion,  find  the  corresponding  point- 
mass  set.  In  practice,  this  means:  given  a finite  set  of  coefficients 
{Anm,  Bnm},  invert  (8)  to  find  the  corresponding  point-mass  parameters 
{Mi,  ri,  <(>i,  Xi>.  The  nature  of  exact  and  approximate  solutions  to 
this  problem  will  be  discussed  in  the  next  section. 


NATURE  OF  SOLUTIONS 


It  is  clear  that  any  practical  point-mass  set  is  defined  by  a 
finite,  bounded  set  of  real  numbers  {Mi,  ri,  <f>i , Xj.}.  Since  this 
bounded  set  gives  rise,  via  (8) , to  an  infinite  number  of  non-zero 
coefficients  {Anm,  Bnmb  the  spherical  harmonic  expansion  corresponding 
to  any  practical  point-mass  set  will,  in  general,  be  infinite.  Tris 
leads  to  a basic  problem  in  point-mass  modeling,  a problem  whose  origins 
will  now  be  elucidated. 


The  discreteness  of  geophysical  data  imposes  a practical  limit* 
on  our  knowledge  of  the  fi>  a structure  of  the  global  gravity  field. 

This  inherent  limitation  prohibits  (long  before  any  computational  or 
time  constraints  would)  the  construction  of  spherical  harmonic  expan- 
sions beyond  some  finite  maximum  degree  and  order.  Thus,  by  attempting 
to  construct  a practical  point-mass  set  from  geophysical  data,  we  are, 
in  effect,  trying  to  match  an  infinite  orthogonal  expansion  with  a 
finite  one.  This  is  a fundamental  difficulty  since  an  exact  match  is 
impossible  (because  all  the  coefficients  cannot  possibly  match) . 

In  the  foregoing,  the  term  "practical"  has  been  emphasized  in  dis- 
cussing point-mass  sets,  since  it  is  only  those  point-mass  sets  defined 
by  bounded  parameters  which  can  be  used  for  practical  purposes.  In 
addition,  there  also  exist  "impractical"  point-mass  sets  (which  will 
be  presently  defined  as  the  limit  of  a sequence  of  practical  sets) . 
These  impractical  sets  are  conceptually  important  as  they  provide  exact 
(though  not  physically  realizable)  solutions  to  the  problem  of  point- 
mass  modeling;  in  comparison,  practical  point-mass  sets  give  physically 
realizable,  though  approximate,  solutions.  Both  cases  will  how  be 
considered  in  more  detail. 


EXACT  SOLUTIONS  (IMPRACTICAL  SETS) 


From  (8) , the  spherical  harmonic  expansion  coefficients  associated 
with  a point -mass  set  are  given  by 


(Anm 

ii 

cos  mXj\ 

I 

= Nnm2^  vi  pi  pnm  (sin  <f>i)  . 

lBnm 

i 

sin  mXjJ 

where 


* In  this  regard  see  Reference  7. 


Consider  the  process  of  taking  the  limits  pi  -*■  0,  Mi  00  such  that 


Mi  P = CV  (10) 

1 l l 

where  the  are  finite  constants;  assume  that  all  the  Pi  are  taken  to 
zero  at  the  same  time,  i.e.  , p^  = 6 and  6 -*■  0.  Then  (10)  can  be 
written 

Mi  = cj  6"n  (ID 

The  Mi  then  diverge  for  any  continuous  or  discrete  sequence  whereby  6 -*■  0. 
For  a particular  value  of  n,  say  n=k,  (11)  becomes 


k .-k 

Mi  = ^ 

Putting  this  into  (9) (using  Pi  = <$)  gives 


~ Nnm  ^ 


y:  C*  Pnm  (sin  <f>i) 


[cos  mXi| 


l sin  mXi f 


fyotl 

- pkm  (sin  <(>i) 

cos  mXi 

Bkm 

i 

sin  mXi 

3.  n > k,  n - k = > 0 


(14) 


(15) 


(16) 


In  the  limit  6 -*  0,  equations  (14),  (15),  and  (16)  become: 

1.  n < k 


= NnmE  C.  Pnm  (sin  *i) 

i l sin  mXf 


2.  n = k 


Akir.  k cos  mXi 

= Nkm^  Ci  Pkm  (sin 
Bkm  i sin  mAi 


3.  n > k 


Up  to  and  including  degree  k,  there  are  (k+1) 2 coefficients,  thus 

(17)  and  (18)  comprise  (k+1) 2 equations.  Since  (10)  defines  the  re- 
lationship between  and  Pi,  there  are  three  independent  variables  for 
each  point-mass.  Thus,  if  I is  the  total  number  of  masses,  then  there 
are  31  independent  variables  in  the  point-mass  set.  Therefore,  the 
minimum  number  of  masses  needed  in  order  that  the  equations  (17)  and 

(18)  be  solvable  is  the  minimum  number  which  contains  at  least  (k+1)*1 
independent  variables.  This  number  is  the  integer  I,  where 


(k+1)2 

< I < 

3 ~ 


(k+1)2 


For  specific  values  of  k,  I is  given  in  Table  1. 

Table  1.  Number  of  Masses  (I)  for  a Particular  Degree  (k) 


APPROXIMATE  SOLUTIONS  (PRACTICAL  SETS) 


If  we  wish  to  use  point-mass  sets  to  approximate  a finite 
spherical  harmonic  expansion,  then  the  radii  of  the  point-masses 
must  be  non-zero;  i.e.,  the  limiting  procedure  p*  = 6 -*■  0 is  only  a 
heuristic  theoretical  device  for  understanding  multipoles.  We  can, 
however,  create  point-mass  sets  whose  corresponding  coefficients 
[via  (9) ] exactly  match  those  of  the  finite  spherical  harmonic 
expansion  (FSHE)  for  degree  n < N (where  N = maximum  degree  of  the 
FSHE)  and  are  sufficiently  small  for  n > N.  This  can  be  done  in  the 
following  manner. 

First,  (17)  and  (18)  can  be  solved  for  each  degree  k <_  N to 

produce  a point-mass  set  {c.  , 4> . , X.}.  Since  C.  = y.  6 , a 

1 1k  1 k -k  1 1 

6 > 0 can  be  chosen,  and  using  y^  = C^  6 , the  corresponding  mass 

can  be  found.  Now,  although  6 > 0,  equations  (17)  and  (18)  are 
still  satisfied;  the  A^  and  for  n > k are  no  longer  identically 
zero,  however,  but  instead 


n > k,  | 

Anm  | 

i | 

rn-k  „ ' n_  , . ,k 

[ - ® Nnm  Z,^iPnm^sln  j 

[cos  mX  . ) 
,£ 

[Bnm  J 

! i 

[sin  mX^j 

By  choosing  6 properly,  the  coefficients  in  equation  (20)  can  be 

made  sufficiently  small;  the  only  limit  comes  from  attempting  to  use 

the  point-mass  sets  so  defined  in  numerical  simulations  because  of 

round-off  error,  etc.  Let  the  coefficients  defined  by  equation  (20) 

k k 

for  each  degree  k be  written  A , B , and  let  them  be  called  the 

nm  nm 

residual  coefficients  of  degree  k. 

Now,  a procedure  for  constructing  point-mass  sets,  whose  intrinsic 

coefficients  for  n N exactly  match  the  given  FSHE  coefficients  and 

for  n > N are  sufficiently  small,  can  be  outlined.  Starting  with 

k = 0 solve  equations  (17)  and  (18) ; using  the  point-mass  set  of  the 

solution,  use  equation  (20)  to  generate  the  A0  , B°  . Subtract 

nm  nm 

these  from  the  remaining  FSHE  coefficients  (0  < n <_  N) , AM,  Bnm,  and 
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now  use  these  subtracted  coefficients  to  solve  for  the  k=l  point-mass 


set,  generating  A'  , B'  from  the  solution.  Then,  subtract  these  from 
nm  nm 

the  remaining  (already  once-subtracted)  coefficients  and  use  the  twice- 
subtracted  coefficients  to  solve  for  the  point-mass  set  corresponding 
to  k=2 , then  k=3,...then  k=N. 

The  number  of  point-masses  needed  to  describe  an  FSHE  of  maximum 
degree  N is  then  the  sum  of  the  number  needed  for  each  degree  k <_ N; 
using  Table  1,  this  number  (call  it  S)  can  be  determined.  The  result 
is  given  in  Table  2. 

Table  2.  Total  Number  of  Masses  (S)  Needed  for 
Expansion  up  to  Degree  and  Order  N 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 

3 

6 

12 

21 

33 

50 

72 

99 

133 

174 

This  result  is  exact  in  the  limit  that  all  6^  0.  (6^  means  that  for 

each  k,  <5  is  chosen  independently  of  the  other  k's.) 

It  may  be  possible  to  use  another  method  for  finding  a point-mass 
set  which  adequately  represents  a given  FSHE,  and,  at  the  same  time, 
has  less  masses  than  the  corresponding  sets  given  in  Table  2.  One  way 
is  to  use  (9)  directly;  in  this  method,  the  p^  would  have  to  be 
small  also,  but  could  be  picked  independently  of  one  another.  Again, 
there  would  be  3 degrees  of  freedom  for  each  point-mass  {p^,  <i>if  . 
The  pj_ , though  independent  of  one  another,  are  not  independent  in  the 
general  sense  because  they  are  constrained  to  be  small  (i.e.,  picked 
beforehand) . 

Thus,  if  the  set  of  non-linear  equations  (9)  is  to  be  solvable, 
there  must  be  enough  variables  {pj_,  <t>i,  A^}  in  the  point-mass  set;  that 
is,  31  (I:  number  of  point-masses)  must  be  greater  than  the  number  of 
coefficients  on  the  left-hand  side  of  (9) . This  number  has  already 
been  determined,  however,  in  Table  1.  Therefore,  it  may  be  possible  to 
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use  41  masses  to  adequately  describe  a (10x10)  FSHE  to  within  an  ac- 
cept«.ble  error;  the  minimum  amount  of  error  is  determined  by  numerical 
limitations  (i.e.,  interaction  between  smallness  of  the  and  a 
particular  computer's  round-off  error,  etc.  One  thing  that  happens  as 
the  -*•  0 in  (9)  is  that  the  Pi  ■>  The  combination  of  small  Pi  and 
large  p^  may  tend  to  further  increase  the  computational  inaccuracy) . 

Since  any  practical  point-mass  set  will,  by  necessity,  only  pro- 
vide an  approximate  model  of  the  geopotential,  it  would  be  useful  to 
have  a method  for  determining  how  "good"  the  approximation  really  is. 
An  appropriate  method  is  defined  and  developed  in  the  next  section. 

GLOBAL  ERROR 


The  geopotential  may  be  represented,  as  has  already  been  remarked, 
by  a finite  spherical  harmonic  expansion 


(sin  c|>) 


cos  mX  + bnm  sin  mX 


] 


(21) 


where  the  a^  and  bju,,  are  derived  from  geopotential  data,  and  N is  the 
maximum  degree  and  order  of  the  expansion. 

The  difference  between  the  geopotential  given  by  (21)  and  the 
point-mass  potential  given  by  (7)  is 

6t  = TG  - T 

_ . 00  n , . n+1  r 

= -5-  £ ^(7)  Pnm  (sin  I (anm  “ Anm)  cos  mX 

n=0  m=°  +(bnm  - Bn*)  sin  mx]  (22) 

where  the  a and  b are  all  zero  for  n > N. 
run  nm 
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GLOBAL  POTENTIAL  ERROR 


b 


The  "global  error"  in  the  potential  will  be  called  ET  and  will  be 


defined  by 

i?2  _ JL 


III 


6T  M dv 


1 

V 


I 


S (r)  r2dr 


(23) 


where 


II 


S (r)  = II  I 6t I 2 da 


(24) 


is  the  integral  of  |6t|  over  a spherical  surface  at  radius  r.  The 


integration  volume  V in  (23)  is  the  volume  in  which  the  point-mass 
method  is  to  be  used;  this  volume  will  be  taken  to  be  the  space  bounded 
by  two  concentric  spheres,  the  outer  one  of  radius  r=b  and  the  inner 
one  of  radius  r=a,  where  b > a >_  R.  Thus,  the  global  error  in  the 
potential  is  seen  to  be  the  root-mean- square  value  of  (22)  in  the 
aforementioned  volume  V,  where 

V = 4ir(b3  - a3) /3 


(25) 


Now,  placing  (22)  into  (24)  and  using  the  orthogonality  properties 
of  spherical  harmonics  (see  Appendix  A)  gives 


■«  ■ o2  So  © 


2n+2 


n n-1 

y'  Nnm 


m=0 


2n+l 


(anm  - ^nm)  + (bnm  “ ®nm) ' 


“'nm 


where  N is  defined  in  (6) , and  where  the  a and  b 

nm  nm  nm 


(26) 
are  zero 


for  n > N. 


Substituting  the  expressions  for  A and  B given  in  (8) 

nm  nm 
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■dMlik., m«I 


into  (26)  and  using  the  addition  theorem  (5)  yields 


S<r>  - c £ (7)  E fSr  + - 2 (anmAnm  + *WW  ] 

n=0  \ / m=0  L J 

4.  rV  V MiMl  V / rir j\n  Pn(cos  ipij) 

*1  r V »-oU  <27, 


Pn (cos  iHj) 
2n+l 


where 


C = 4tt 


and  where  is  the  angle  between  the  position  vectors  of  the  i*-h  a 
jth  point-masses. 

Expression  (27)  for  the  surface  integral  S(r)  can  be  put  into 
closed-form  by  utilizing  the  first  result  of  Appendix  B;  thus 

S (r)  = C X.  2n+l  f3^10  + bnm  ~ 2 ^ anmAnm  + bnmBnm)J 

r»  = n ' ' m=D  J 


* Z ^ p(ctij'  6lj>  <29) 

1 3 r^Vr^r  j 

where  F(a,  6)  is  the  elliptic  integral  of  the  first  kind  (see  Appen- 
dix C)  with 


= arccos 


r - r^r-i 


rz  + ri] 


Bil  = cos  (*5^ii) 


The  global  potential  error  can  now  be  evaluated.  Recalling 


E2  = 


f b 

I S (r)  r^ 


and  placing  (27)  into  it  [ (27)  is  more  useful  in  the  present  context 
than  (29) ] yields 


n -1 


T v ^ I>i  (a,b,n)  £ f^nm  + bnm ( anm^W.  + bnmBnmj] 

n=0  m=0  L -> 

, CR3\p  tp  MiMj  T i v.  >> 

+ — 44-  — T~  b2  (a,b,n) 


M 


where 


1 T/R\2n*l  /R  \2  n- 1~| 

Ll(a,b,n)  j 

-b.nl  = T Pn  (cos  *i-j) 

n=0  R [\  a / V b / J (2n-l)  ("2n+i 


) (33) 


Since 


(2n-l) 


1 = —L — i L_\ 

(2n+l)  2 \2n-l  2n+l) 


(33)  can  be  written 


L2  (a,b,n)  = -li 


rrr  Ki|  — > ^ii)  - k2/  * 


2R  rir 


where  the  functions  Kj  and  K2  are  given  in  Appendix  B,  and 

f(r)  |*  = f(y)  - f(x)  (35) 

Equation  (31)  is  a rather  complicated  expression  for  the  global 
potential  error,  primarily  because  of  the  function  L2.  In  practice, 
however,  (31)  will  be  of  minimal  importance  since  the  usual  aim  of 
poinc-mass  modeling  is  to  match  the  gradient  of  the  geopotential 
rather  than  the  potential  itself.  To  this  end,  a "global  gradient 
error"  will  now  be  derived. 


GLOBAL  GRADIENT  ERROR 


The  global  gradient  error  E will  be  defined  as  follows 

G 


EG  = ? fff  l6VTl2  dv 


where  the  difference  in  the  gradients  of  the  point-mass  potential 
(7)  and  the  geopotential  (21)  is 

6VT  = VT  - VT  = V(T  - T)  = V6T  (2 

G G 

Thus,  (36)  can  be  written 


■ s III' 


V6T  • V6T  dV 


[Here,  the  volume  V is  the  same  as  was  used  in  defining  the  global  po- 
tential error  (23) • ] 

Equation  (38)  can  be  put  in  a more  useful  form  by  utilizing 
Green's  First  Identity  (Reference  6,  p.  11) : if  u and  W are  any  two 
scalar  functions  (potentials) , then 


fff (U?2w 


+ VU  • VW) dV 


= ffu  n* 


where  S is  the  closed  surface  bounding  V and  n is  a unit  vector  normal 
to  the  surface  S and  directed  out  of  the  volume  V.  Setting  U = W = 6T 
and  using  V26t  = 0 (since  T and  T in  (22)  are  both  harmonic  in  V) 

G 

allows  (38)  to  be  written 

EG  = V **  * dS  (40) 

S 

The  surface  S in  (40)  consists  of  the  inner  and  outer  spheres  which 

bound  V;  thus,  the  normal  on  the  outer  sphere  (r  = b)  is 

nfa  = r = r/|r|,  while  on  the  inner  sphere  (r  = a)  the  normal  is 

A 

n = -r. 
a 


Therefore, 


% ’ v = -na 
and  (40)  becomes 


E2  = — 
G V 


fl 

cb 


■ if- 

hff  >6T|Z  da 

cr 


V = r • V = 

3r 


6T  ( -^  6t  1 b2do 


(41) 


If 


<5T  ( — 6t  ) a^do 


Z-  6T  rzda 


(42) 


Finally,  using  (24) , the  square  of  the  global  gradient  error  becomes 


eg  = lrs(r) 


where  S (r)  is  given  by  (29). 

Evaluating  the  partial  derivatives  in  equation  (43)  gives 


(43) 


r 


i 


and  where 


Hij (a,b)  = 


COS  j 


b 

(46) 


Here  F(a,8)  is  again  the  elliptic  integral  of  the  first  kind  (Appen- 
dix C)  with 


= arccos  / r rir  j 

13  \ y-2j 


+rir  j 


(47) 


6ii  = cos  (‘jiJ’i-j) 


■'ij  _ v J*io 
and  where,  as  in  (42) , use  has  been  made  of  (35) . 

The  analytic  expression  (44)  for  the  square  of  the  global  gradient 
error  can  be  written  in  a more  compact  form  by  again  using  the  follow- 
ing change  of  parameters 


Mi 

M 


pi 


(48) 


(the  latitude  and  longitude  parameters  <(>^  and  Ai  remain  unchanged) 
Using  (48) , (44)  becomes 

EG  = c>  (*  ~ ^i^i  + J pipjhij\ 


1 3 


where 


C’  = 


CR  3G2R3 


(49) 


(50) 


b3-a3 


(from  (25),  (28)  and  G = kM/R2,  where  G:  mean  gravity  on  the  earth's 
surface)  and  where 


f = Y G (a,b)  X N-1  (a2  + b2  ) 

^ n “ nm\  nm  nm/ 

n=0  — 


(51) 


m=0 
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g.  = 2 Y.  Gn(a,b)p”  £ pnm(sin  <f>i>  (anm  cos  mXi  + bnm  sin  mXi j 
1 n=0  m=0  V 


(52) 


(PiPj)  F(aij'  B 


. .) 


13 


(pVfp.p.)2  - 2p2p1p_.  cos 


a/R 

|p=b/R 


(53) 


with 


/p2-PiP-i 

= arccos  i x J 


(P2+PiPj 


(54) 


3i:j  = cos 


iy 


(55) 


Here  G (a,b)  is  given  by  (45)  and  F(a,6)  is  again  an  elliptic 
n 


integral  of  the  first  kind. 


OPTIMIZATION  PROCEDURE 


Now  that  an  analytic  expression  (49)  has  been  obtained  for  the 
global  gradient  error,  it  will  be  possible  to  develop  an  optimization 
procedure  for  point-mass  sets.  The  obvious  requirement  for  optimiza- 
tion is  that  the  global  gradient  error.  Eg,  be  a minimum  with  respect 
to  the  independent  parameters  of  the  point-mass  set.  Since  E^  is 
always  positive,  its  minima  occur  for  the  same  values  as  those  of 


E2;  thus,  the  optimum  point-mass  parameters  are  to  be  found  by  solving 
G 


the  following  equations; 


3 E i 


3 P, 


= 0 


7,  E2  = 0 


k G 


(56) 
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r — . 


where  the  gradient  with  respect  to  the  coordinates  of  the  k*"*1  point- 
mass  is 

„ 3 , 1_  3 ' 1 3 

Vk  = rk  3rk  + k rk  3^  + k rR  cos  9*k 


1 / - 9 , x 1 8 

= - | r +4>.  — rr—  + X 


r l k 3pk  k Pk  ' k Pk  cos  0>k  3X1 


(57) 


where  rk,  and  X^  are  unit  vectors  in  the  direction  of  greatest 

increase  of  r^,  (J)^,  and  Xj,,  respectively. 

-.2 


Applying  the  conditions  (56)  to  the  expression  (49)  for 


gives 


gk  ^yjhkj 
J 

V g = Vy.Vh  . 
k k V ] k k] 


(58) 


(59) 


(where  it  must  be  understood  that  Vk  h^  = lim(j  -*•  k)  Vk  h^  and  not 
Vk  hkk  = ^k  kjh^).  Assuming  that  {h^}  is  a non-singular 

matrix  allows  equation  (58)  to  be  inverted  to  give 


j . = Y'  h . * g. 

1 Y ]k  k 


(60) 


Placing  (60)  into  (59)  gives  an  optimizing  equation  which  does  not 
contain  any  explicit  masses  1% 


7kgk  = hji  [Vkhkj]gi 


(61) 


[Remember  that  (61)  represents  three  equations  for  each  mass,  because 
of  (57) . Expressions  for  V h . are  given  in  Appendix  E. ] 

The  derivation  of  the  optimizing  equations  (58)  through  (61)  is 
much  more  straight-forward  than  their  solution.  If  equation  (61) 
could  be  solved  exactly  to  give  the  optimal  positions,  then  (60) 
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could  be  easily  used  to  complete  the  set  by  producing  the  optimum 
mass  values.  (61) , however,  is  highly  non-linear  and  a solution  can 
probably  be  obtained  only  by  an  extensive  numerical  search  algorithm; 
the  speed  of  computation  and  accuracy  of  a final  solution  are  diffi- 
cult to  estimate,  primarily  because  of  the  presence  of  h^!  in 
(61). 

Another  approach  is  to  try  to  find  the  minimum  of  E2  numerically, 

G 

with  the  aid  of  the  optimizing  equations,  rather  than  trying  to 
solve  (61)  directly  (either  numerically  or  analytically) . In  parti- 
cular, for  an  arbitrary  set  of  positions,  the  h^_.  can  be  calculated 
and  (assuming  the  matrix  {h^}  is  non-singular)  (60)  can  be  readily 
utilized  to  compute  the  optimum  choice  of  masses  for  the  arbitrary 

positions.  (That  these  masses  yield  a minimum  for  E2  follows  from 

G 

o 

the  fact  that  is  positive  and  quadratic  in  the  u^.  Let  the 
masses  which  satisfy  (60)  for  arbitrary  locations  be  denoted  by  P^; 

o , 

then  can  be  written 

EG(V  = EG(V  + Q (62) 

where  Q is  the  quadratic  form 


Q = ? f hijEi£j 


and  y = y , + £..  Since  E2  > 0,  then  Q > - E2(p  ) ; and  since  Q can 
K K 1 G G K 

always  be  written  in  canonical  form®,  i.e., 

(64) 


then  if  any  ai  < 0,  there  exists  a direction  in  which  Q will  grow 

unboundedly  large  and  negative  exceeding  the  lower  bound  - E2  (p  ) . 

G R 

Therefore,  all  a^  must  be  non- negative  from  which  it  follows  that 
Q >_  0,  and  thus  the  specify  the  minimum.) 
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Once  the  p have  been  found,  via  (60) , for  a given  set  of  posi- 
tions, a gradient  descent  method  can  be  utilized  to  seek  out  the 

minimum  of  E?,  with  respect  to  the  position  parameters.  (Since  the 

G 

q 

gradient  descent  method  is  thoroughly  discussed  elsewhere-,  it  will 
not  be  considered  in  detail  here. ) The  standard  gradient  descent 
algorithms  need  only  be  appended  in  a simple  manner:  at  the  end  of 
each  iteration  use  (60)  to  determine  optimum  masses  for  the  new 
(iterated)  positions:  then  the  descent  algorithm  need  work  on  31 
variables  instead  of  41,  where  I is  the  number  of  masses. 

It  is  important  to  note  that  care  must  be  taken  in  choosing  the 
number  of  masses  for  use  with  this  optimization  procedure.  As  was 
shown  earlier,  each  finite  spherical  harmonic  expansion  of  degree  and 
order  N has  associated  with  it  an  exact  solution  (as  a limiting 
case) , and  each  exact  solution  has  its  own  particular  number  of 
masses  S.  Therefore,  if  the  degree  and  order  of  the  FSHE , which  is 
being  approximated,  is  N and  the  point-mass  set,  which  initializes 
the  numerical  optimization  procedure,  consists  of  S or  more  masses, 
the  search  algorithm  may  try  to  reach  the  limiting  values  r^  -*■  0 and 
->  °o  (i  =1,...,S);  this  will  clearly  cause  problems. 

In  practice,  however,  it  is  often  useful  to  have  as  few  masses 
as  necessary.  In  particular,  it  has  been  determined  at  NAVSWC 
that  the  time  required  to  compute  a gravity  vector  (on  the  Mk  98  Mod  0 
fire  control  computer)  is  the  same  for  a 9x9  FSHE  as  for  a point-mass 
expansion  consisting  of  about  50  masses.  In  this  case,  it  is  obvious 
that  fewer  than  50  masses  are  needed  if  computation  time  is  to  be 
conserved  (an  important  consideration  in  fire  control) . Thus,  the 
availability  of  an  efficient  point-mass  optimization  procedure  can 
greatly  enhance  existing  fire  control  capabilities. 
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DISCUSSION 


In  order  to  give  some  credence  to  the  theoretical  claims  contained 
in  this  report,  the  following  numerical  wo^-k  was  done.  Since  equation 
(8)  is  linear  in  the  mass  values  M^,  it  can  be  inverted  to  yield  these 
mass  values  provided  certain  conditions  are  met:  1)  there  are  as  many 
masses  as  coefficients  Am,,,  Bnm,  2)  the  positions  (ri,  4>i,  Aj_)  are 
chosen  beforehand  and  in  such  a manner  that  the  matrix  which  must  be 
inverted  is  non- singular . Also,  in  order  to  achieve  a form  appropriate 
for  inversion,  the  double  indices  n,m  appearing  on  both  sides  of 
equation  (8)  had  to  be  "stretched"  into  a single  index  (suitable  for 
describing  a column  vector) . 

The  coefficients  used  for  A^,,  Bnm  came  from  the  NASA  GEM-6  set 
(n=0,l  terms  equal  to  zero)  with  maximum  n,m  equal  to  nine;  thus, 
there  were  100  coefficients.  The  angular  positions  of  the  required 
100  masses  were  chosen  at  random  and  the  radii  of  these  point-mass 
positions  was  the  same  for  all  masses  (except  for  one  which  was  always 
placed  at  r=0) . The  value  of  the  common  radius  was  set  at  .15,  .30, 
.45,  .60,  .75,  and  .90  earth  radii  (R) . For  each  of  these  values, 
equation  (8)  was  inverted  and  the  corresponding  mass  values  were  found 
(the  angular  positions  were  the  same  for  all  radii) . 


To  check  how  well  each  of  these  point-mass  sets  approximated  the 
original  finite  spherical  harmonic  expansion  (9x9) , both  expansions 
(point-mass  and  FSHE)  were  used  to  compute  the  values  of  the  gravity 
disturbance  vector  at  100  randomly  chosen  positions  between  1.0  and  1.01 
earth  radii.  The  root-mean- square  (RMS)  differences  were  then  calcu- 
lated for  each  of  the  point-mass  set  radii  .15,...,  .90,  and  the 
results  are  shown  in  Table  3. 


22 


I 


Table  3.  RMS  Difference  in  Components  of 
Gravity  Disturbance  Vector  T 


RMS  of  6 Components  (Milligals) 


Radius  (R) 

6r 

6<}i 

6X 

.15 

4.22 

2.68 

14.3 

.30 

18.9 

13.7 

19.2 

.45 

64.0 

42.2 

37.4 

.60 

165. 

99.1 

106. 

.75 

1060 

774. 

905. 

.90 

40400 

29700 

24500 

Although  the  analytic  expression  for  the  global  gradient  error  was 
not  used,  these  results  still  provide  quantitative  verification  of  some 
of  the  theoretical  developments  of  this  report.  In  particular,  they 
support  the  contention  contained  in  the  section  on  "Nature  of  Solutions", 
that  the  point-mass  sets  provide  for  better  approximations  the  closer 
the  masses  are  to  the  center  of  the  earth.  The  optimal  point-masses 
should,  however,  have  radial  values  which  do  not  tend  completely  to  zero 
(though  small)  values.  (Otherwise,  these  optimal  point-mass  sets  would 
not  exist,  in  a practical  sense,  as  has  already  been  discussed.) 


CONCLUSION 


In  this  report,  the  basic  theory  of  global  point-mass  modeling 
has  been  discussed.  Since  all  the  information  obtainable  from  discrete 
global  geophysical  data  can  be  essentially  incorporated  into  the  coef- 
ficients of  a finite  spherical  harmonic  expansion,  the  independent 
parameters  of  a point-mass  expansion  are  effectively  functions  of  these 
coefficients.  The  implicit  functional  form  was  determined  and  is  given 
by  (8) . This  functional  relationship  showed  that  a practical  point- 
mass  expansion  will  have  an  associated  infinite  spherical  harmonic 
expansion.  This,  in  turn,  exposed  a basic  problem  in  point-mass  modeling, 


i 
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that  of  matching  an  infinite  with  a finite  orthogonal  expansion.  The 
functional  relationship  then  allowed  the  nature  of  solutions  (both 
exact  and  approximate)  for  this  basic  problem  to  be  carefully  considered. 

In  order  to  better  gauge  how  well  a particular  point-mass  set 
approximates  the  global  geophysical  data,  the  concept  of  "global  error" 
was  defined  and  developed.  A desire  to  minimize  the  global  gradient 
error  led,  in  turn,  to  a set  of  simultaneous  non-linear  equations  whose 
solution  allows  for  the  determination  of  the  optimal  positions  and 
values  of  a given  number  of  point-masses. 

Although  a procedure  for  solving  these  equations  was  sketched,  it 
is  well  recognized  "that  the  solution  of  simultaneous  non-linear  equa- 
tions is  usually  a very  difficult  problem"10,  a problem  beyond  the  scope 
of  the  present  work.  Thus,  the  implementation  of  the  methods  described 
herein  are  left  as  a (hopefully)  well-posed  numerical  problem,  content- 
ing ourselves  at  present  with  their  derivation. 


ORTHOGONALITY  PROPERTIES  OF  SPHERICAL  HARMONICS 


The  following  results  will  be  written  in  the  notation  of  Heiskanen 
and  Moritz  (Ref.  6,  p.  29). 


1.  Spherical  Harmonics 


R (0,  1)  = P (cos  0)  cos  mX 
nm  nm 


S (0,  X)  = P (cos  0)  sin  mX 
nm  nm 


(65) 

(66) 


2.  Orthogonality  Relations 


If 


R (0,  X)  R (6,  X)  do  = e yn-+mj  ’■  6 <5 

nm  rs  2n+l  m (n-m)  ! nr  ms 


(67) 


If 


S (0,  X)  s (0,  X)  do  = e ' 6 6 

nm  rs  2n+l  m (n-m) ! nr  ms 


(68) 


JJ  R^m  (0,  X)  Srs  (0,  X)  do  = 0 , always 


(69) 


Here,  6nr  is  the  Kronecker  delta, 


1 , n = r 


Jnr 


0 , n f r (70) 

and  em  = *5(1  + 6mo)  . 

The  coordinate  system  employed  in  the  above  definitions  and 
results  is  a spherical  polar  coordinate  system: 

*2ir  /-if 


iu  - r s 

— \ f\ 


sin  0d0dX 


(71) 


X=0  0=0 


A-l 


The  relation  of  spherical  polar  coordinates  6 and  X and  the 
geocentric  coordinates  <f>  and  X is 


" 6 " 2 


(72) 


with  the  azimuthal  angle  X being  the  same  in  both  systems.  Thus, 
cos  6 = sin  4>  (73) 

and 


•2tt  /»tt_ 
2 


IhLf. 


cos  <f>di>dX 


(74) 


x=o  - 


APPENDIX  B 

CLOSED  FORM  EXPRESSIONS  FOR  CERTAIN  INFINITE 
SERIES  OF  LEGENDRE  POLYNOMIALS 
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CLOSED  FORM  EXPRESSIONS  FOR  CERTAIN  INFINITE 
SERIES  OF  LEGENDRE  POLYNOMIALS 


Result: 


Kj  (z , ip) 


CO 


E 

n=0 


1_ 

z2n+l 


Pn  (cos  Ijj) 
2n+l 


*jF  (a,  6) 


(75) 


where  F(a,  6)  is  the  elliptic  integral  of  the  first  kind  (see  Appen- 
dix C)  with 

/ z2-l 

a = arccos  

\z2+l 

6 = cos(<|j/2) 

Derivation: 


x~2n+2'  ' pn(cos  M dx 
dx 

2 cos  ip  + x4)*5  (76) 

where  (3)  and  (4)  have  been  used. 

Utilizing  formula  3.165-2  of  Reference  11  leads  to  (75). 


Result : 


k2(z,  IP)  = £ 

n=0  z2n+1  2n  - 1 


= -1  + z£l  - z2(1-2z2  cos  \()  + z4) 

+2^cos  ip  J4  (z,  1,  - cos  i/j)  - J2  (z,  1, 


- cos  ||/)  (77) 


B-l 


where  J2  and  J4  are  given  by  (97)  and  (99) , respectively. 


Derivation: 


K2(z,  i|»)  =£ 


1 Pn (cos  ^) 
,2n-l  2n  - 1 


/•oo 

I 00 

* E 

J z n=l  x2n 


(cos  \|>)  dx 


ft 


= -1+1  x2  (1  - 2xz  cos 


+ X4)"15  - l] 


dx  (78) 


where  the  last  step  follows  from  (3)  and  (4) . 

Integrating  (78)  by  parts  gives 
K2(z,  i|/)  = -1  + z ("l  - z2(l  - 2z2  cos  iji  + z4) 


K2(z,  \p)  = -1  + z £l  - z2(l  - 2z2  cos  iji  + z4) 

/*°°  4 , 2 

+ f x4.,cosJ>  dx 

(1  - 2x2  cos  Ip  + x4)*2 


+ x4)3^ 


(1  - 2x*  cos  <p  + x ) 


using  (97)  and  (99)  of  Appendix  D then  gives  the  result  (77) 
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THE  FIRST  AND  SECOND  KINDS 
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SOME  PROPERTIES  OF  ELLIPTIC  INTEGRALS  OF 
THE  FIRST  AND  SECOND  KINDS 


l 


■ 


The  following  results  come  from  Section  8.1  of  Reference  11. 
1.  Elliptic  integral  of  the  first  kind 


F(4>,  k)  = 


/' 

^ o 


da 


(80) 


Vl-k2  sin^  a 


2.  Elliptic  integral  of  the  second  kind 

E (<f> , k)  = 


r4> 

- IT- 

J Q 


k2  sin2  a da 


(81) 


3.  Derivatives 


9F 

3$ 


^ 1-k 


? • 2 

sm 


(82) 


9 4>  v 


1-k2  sin2 


3F 

3k 


1 E 


k,2F 


k sin  d>  cos 


(83) 

(84) 


k2  sin2 


3E  _ E - F 
3k  k 


(85) 


where 


k'  =\l-) 


Elliptic  integrals  of  the  first  and  second  kinds  also  have  the 


12 


following  special  values. 

F ( <}> , o)  = E (<(> , o)  = 


(86) 
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APPENDIX  D 

EVALUATION  OF  CERTAIN  INTEGRALS  WHOSE  INTEGRANDS 
CONTAIN  THE  FACTOR  (x4+2b2x2+a2 ) 


EVALUATION  OF  CERTAIN  INTEGRALS  WHOSE  INTEGRANDS 
CONTAIN  THE  FACTOR  (x4+2b2x2+a2 ) "^ 


Formula  3.165-2  of  Reference  11  is 


T.x(u,aS  b^) 


2 = 


/•oo 


dx 


4+2b2x2+a4 


= F (a , r) 
2a 


(89) 


\a  > b^  > az  > 


0,  U > o] 


where  F(a,  r)  is  an  elliptic  integral  of  the  first  kind  (see  Appen- 
dix B)  and 

(90) 


a = arccos 


( u2  - a2) 
Vu2  + a2} 


r = 


Va2  - b2 


a V 2 


Consider  the  following: 


J x (u,a2  , b 2)=^ 


dx 


(x4+2b2x2+a4)  ^ 


1 3lx 
2a3  3a 


“T  r)  " TT  7“  F(a'  r) 
4a 3 4a4  8a 


4a4 


1 , . 3F  8a 

— F (a,  r)  - — — 
a 3a  3a 


3a  _ 3F  3r"j 
3a  3r  3a  I 


(91) 


(92) 


Differentiating  (90)  and  (91)  with  respect  to  a yields 


3a  2u 


(93) 


3a  u2+a2 


da 


VTa2 


Now,  using  (93)  and  (94)  along  with  (82)  and  (84)  of  Appendix  C 
(92)  becomes 


Jj (u,a2 ,b2)  = — F(a,r) 
4a3 


4aH 


2u 


VT: 

where  r'  = ^1-r2 . 


^l-r2  sin2a  (u2+a2) 

b2 |f (a,  r)  - r'2  F(a,  r) 

r’2  a2‘VS^brL  r 


r sin  a cos  a 


■yi-r2sin2a 


(95) 


Using  (90)  and  (91) , (95)  can  also  be  written 


Jl  (u,a2,b2)  = ~r 
4a° 


2b^ 

I^Flc"  r)  '1 7^7  El“'  r) 


2ua (u2+2b2+a2) 


(96) 


(a2+b2)  (u2+a2)  \u‘*+2b2u2+a1< 


Formula  3.165-6  of  Reference  11  gives 


J2 (u,a2 ,b2)  = 


f 


x2dx 


(x4+2b2x2+a2)  ^ 


2 (a4-b4) 


E(a,  r)  - 


4a (a2 -b2) 


F (a,  r) 


u(u2-a2) 


(97) 


2 (a2+b2) (u2+a2) Vu4+2b2u2+a2 
fa2  > b2  > - «,  a2  > 0,  u > ol 


t 


Also,  formula  3.165-7  of  Reference  11  is 
J3 (u,a2,b2)  = 


-f 

•'ll 


(x2-a2)2  dx 


(x4+2b2x2+a4)V;! 


a2-b2 


^F (a,  r)  - E(a,  r)J 


2 2 

u -a^ 


(98) 


u2+a2  Vu4+2b2u2+a4 
b2  | < a2 , u>o] 

Here,  and  in  (97),  a and  r are  still  defined  by  (90)  and  (91), 
respectively. 

Finally,  (96) , (97) , and  (98)  can  be  used  to  write 
x4dx 


Ji^Uja^b2)  = 


f 

-'ll 


(x4+2b2x2+a4) ^ 


= J3  + 2a2 J2  - a4 Jj 


3 a , b2 (a+2b)  , 

F (a,  r)  - — 1 — E(a,  r) 


4 a2-b2 


2 (a4-b4) 


+ u[b2(u2-2a2)  - *sa2(u2+a2)] 
(a2+b2) (u2+a2)  ^i4+2bu2+a4 

[|b2|  < a2,  u > 0] 


(99) 
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DERIVATIVES  OF  THE  FUNCTIONS  hij 


3h 


^3 

3pi 


1 

2Pi 


~ ( P iP j ) F(aij,  Bij) 


P 

+ 


p*1  + 2(pipj)2  - 3p2pipj  cos  i^ij 


P4  + (PiP  j) 2 - 2p2PiPj  cos  ^ij  1 /2 


a/R 


(100) 

p=b/R 


3h^ 
8<f) 


ii  = Dij  |cos  <|>i  sin  <f> j - sin  <£i  cos  4> j ccs(Ai  - Aj)J 


(101) 


3hij 


3^i 


= j j^-cos  4>i  cos  4>  j sin(Ai  - A j ) J 


(102) 


where 


a^j  = arccos 


PjPj 


P + PiP j 


(103) 


Bij  = cos^ij) 


(104) 


cos  = sin  sin  <{> j + cos  <fi  cos  <J> j cos(Ai  - Aj)  (105) 


Di).i 


(PiP j) 

4(1-6": . ) 


ID 


E(cij,  6i j)  — ( 1— 3 i j ) F(cij,  8 i j ) sinaij  cosoiij 


b2 


ID 


\l  1-62  . sin2a . . 

" n 


P3PiPj 


[p4  + (pipj)2  - 2p2pipj  cos 


a/R 

p=b/R  (106) 


Note:  the  relation 
3h 


ii  3hij 

= lim 


(107) 


3pi  j .+  i 9Pi 


E-l 


I 1 
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